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Abstract 

We discuss a new class of driven lattice gas obtained by coupling the one-dimensional totally asymmetric simple exclusion 
process to Langmuir kinetics. In the limit where these dynamics are competing, the resulting non-conserved flow of particles on 
the lattice leads to stationary regimes for large but finite systems. We observe unexpected properties such as localized boundaries 
(domain walls) that separate coexisting regions of low and high density of particles (phase coexistence) . A rich phase diagram, 
with high an low density phases, two and three phase coexistence regions and a boundary independent "Meissner" phase is 
found. We rationalize the average density and current profiles obtained from simulations within a mean-field approach in the 
continuum limit. The ensuing analytic solution is expressed in terms of Lambert W-functions. It allows to fully describe the 
phase diagram and extract unusual mean-field exponents that characterize critical properties of the domain wall. Based on the 
same approach, we provide an explanation of the localization phenomenon. Finally, we elucidate phenomena that go beyond 
mean-field such as the scaling properties of the domain wall. 

PACS numbers: 02.50.Ey, 05.40.-a, 64.60.-i, 72.70.+m 



I. INTRODUCTION 



Many natural phenomena driven by some external 
field or containing self-propelled particles evolve into sta- 
tionary states carrying a steady current. Such states 
are characterized by a constant gain or loss of energy, 
which distinguishes them from thermal equilibria. Exam- 
ples range from biological systems like ribosomes moving 
along m-RNA or motor molecules "walking" along molec- 
ular tracks to ions diffusing along narrow channels, or 
even cars proceeding on highways. In order to elucidate 
the nature of such non-equilibrium steady states a vari- 
ety of driven lattice gas models have been introduced and 
studied extensively |l| . Here we focus on one-dimensional 
(ID) models, where particles preferentially move in one 
direction. In this context, the Totally Asymmetric Sim- 
ple Exclusion Process (TASEP) has become one of the 
paradigms of non-equilibrium physics (for a review see 
Refs. US Si). In this model a single species of parti- 
cles is hopping unidirectionally and with a uniform rate 
along a ID lattice. The only interaction between the par- 
ticles is hard-core repulsion, which prevents more than 
one particle from occupying the same site on the lattice; 
see Fig. ^ 

It has been found that the nature of the non- 
equilibrium steady state of the TASEP depends sensi- 
tively on the boundary conditions. For periodic bound- 
ary conditions the system reaches a steady state of con- 
stant density. Interestingly, density fluctuations are 
found to spread faster than diffusively 0- This can be 



understood by an exact mapping to a growing in- 
terface model, whose dynamics in the continuum limit 
is described in terms of the KPZ equation Q and its 
cousin, the noisy Burgers equation 9]. In contrast to 
such ring systems, open systems with particle reservoirs 
at the ends exhibit pha se transitions upon varying the 
boundary conditions [lOj • This is genuinely different from 
thermal equilibrium systems where boundary effects usu- 
ally do not affect the bulk behavior and become negligible 
if the system is large enough. In addition, general the- 
orems do not even allow equilibrium phase transitions 
in one-dimensional systems at finite temperatures (if the 
interactions are not too long-range) [Tl| . 
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FIG. 1: Illustration of the Totally Asymmetric Simple Exclu- 
sion Process with open boundaries. The entrance and exit 
rates at the left and right end of the one-dimensional lattice 
are given by a and /3, respectively. 

Yet another difference between equilibrium and non- 
equilibrium processes can be clearly seen on the level 
of its dynamics. If transition rates between microscopic 
configurations are obeying detailed balance the system is 
guaranteed to evolve into thermal equilibrium 61]. Sys- 
tems lacking detailed balance may still reach a steady 
state, but at present there are no universal concepts like 
the Boltzmann-Gibbs ensemble theory for characterizing 
such non-equilibrium steady states. In most instances 
one has to resort to solving nothing less than its full dy- 
namics. It is only recently, that exact (non-local) free 
energy functionals for driven diffusive systems have been 
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derived [HEl. 

This has to be contrasted with dynamic processes such 
as the adsorption-desorption kinetics of particles on a lat- 
tice coupled to a bulk reservoir ( "Langmuir Kinetics" , 
LK), see Fig- EI Here, particles adsorb at an empty site 
or desorb from an occupied one. Microscopic reversibility 
demands that the corresponding kinetic rates obey de- 
tailed balance such that the system evolves into an equi- 
librium steady state, which is well described within stan- 
dard concepts of equilibrium statistical mechanics. If in- 
teractions between the particles other than the hard-core 
repulsion are neglected, the equilibrium density is solely 
determined by the ratio of the two kinetic rates 0], as 
given by the Gibbs ensemble. 
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FIG. 2: Illustration of Langmuir kinetics, lua and ujd denote 
the local attachment and detachment rates. 

The TASEP and LK can be considered as two of the 
simplest paradigms which contrast equilibrium and non- 
equilibrium dynamics and stationary states. Langmuir 
kinetics evolves into a steady state well described in terms 
of standard concepts of equilibrium statistical mechanics. 
Driven lattice gases such as the TASEP evolve into a sta- 
tionary non-equilibrium state carrying a finite conserved 
current. Whereas such non-equilibrium steady states are 
quite sensitive to changes in the boundary conditions, 
equilibrium steady states are very robust to such changes 
and dominated by the bulk dynamics. In the TASEP the 
number of particles is conserved in the bulk of the one- 
dimensional lattice. It is only through the particle reser- 
voirs at the system boundaries that particles can enter or 
leave the system. In LK particle number is not conserved 
in the bulk. Particles can enter or leave the system at 
any site. Depending on whether we consider a canonical 
or grand canonical ensemble the lattice is connected to 
a finite or infinite particle reservoir. Unlike the steady 
state of the TASEP, the equilibrium steady state of LK 
does not have any spatial correlations. 

Combining both of these processes may at first sight 
seem a trivial exercise since one might expect bulk effects 
to be predominant in the thermodynamic limit. This is 
indeed the case for attachment and detachment rates, 
loa and ujd, which are independent of system size N. 
For large but finite systems interesting effects can only 
be expected if the kinetics from the TASEP and LK com- 
pete. Then, as we have shown recently novel behav- 
ior different from both LK and the TASEP appears. In 
particular, one observes phase separation into a high and 
low density domain for an extended region in parameter 
space. 

When should one expect competition between bulk 
dynamics (LK) and boundary induced non-equilibrium 



effects (TASEP)? Let's consider the following heuristic 
argument. A given particle will typically spend a time 
t ~ 1/ujd on the lattice before detaching. During this 
"residence" time the number of sites n explored by the 
particle is of the order of n ~ t. Hence, for fixed ujd, 
the fraction n/N ~ 1/(u>dN) of sites visited by a par- 
ticle during its walk on the lattice would go to zero as 
N — > oo. Only if we introduce a "total" detachment rate 
by Qd = Nlod and keep it constant instead of lod as 
N — > co will the particle travel a finite fraction of the 
total lattice size. Similar arguments show that a vacancy 
visits an extensive number of sites until it is filled by at- 
tachment of a particle if lua scales to zero as loa = CIa/N 
with a fixed "total" attachment rate ft a- In other words, 
competition will be expected only if the particles live long 
enough such that their internal dynamics or the external 
driving force transports them a finite fraction along the 
lattice before detaching. Then, particles spend enough 
time on the lattice to "feel" their mutual interaction and, 
eventually, produce collective effects. In summary, com- 
petition between bulk and boundary dynamics in large 
systems (N 3> 1) is expected if the kinetic rates uja and 
lod decrease with increasing system size N such that the 
total rates £Ia and £Id with 

n A = Lo A N, n D = uj D N (i) 

are kept constant with N. 

The competition between boundary and bulk dynamics 
is a physical process that has, to our knowledge, not yet 
been studied in the context of driven diffusive systems. 
In previous models emphasis was put on the analysis of 
boundary induced phenomena in driv en g ases of mono- or 
multi-species of particles E IHIlrilr^ . ll8lll9| . in presence 
of interactions (see e.g. R ef. | 2fl l2l | ). disorder [22,[2i| or 
local inhomogeneities |24Ll25ll, particles with sizes larger 
than the lattice spacing |26t |27|. lattices with different 
geometries (e.g. multi-lanes lattice gases j^), or systems 
in presence of several conservation laws (for a review see 
e.g. Ref. [3). 

In this manuscript we explore the consequences of par- 
ticle exchange with a reservoir along the track (LK) on 
the stationary density and current profiles and the en- 
suing phase diagram of the TASEP. A short account of 
our ideas has been given recently 0], where we have 
introduced the model and have shown how our Monte 
Carlo results can be rationalized on the basis of a mean- 
field theory, which we also solved analytically. The pur- 
pose of the present manuscript is to give a complete and 
comprehensive discussion of the topic. We will present 
results from Monte-Carlo results for the full parameter 
range of the model including the particular case where 
on- and off-rates equal each other, which were left out in 
our short contribution [l^ due to the lack of space. In 
addition, we will give the full reasoning for the derivation 
and analytical solutions of our mean-field theory. Here, 
additional insight is gained by identifying a branching 
point that explains all the features of the density pro- 
files and phase diagram analytically. In particular, we 



2 



show that a new critical point organizes the topology of 
the diagram and leads to unexpected phenomena already 
briefly discussed in Ref. 0. In a recent work, Evans et 
al. [56j have rephrased the mean-field analysis first given 
in Ref. 0] and reproduced some of our results. The 
mean-field equations are, however, left in their implicit 
form and thus miss the interesting features we will obtain 
from the identification of a branching point. 

These effects differ from those known in reference mod- 
els of equilibrium and non-equilibrium statistical me- 
chanics like LK and the TASEP. Indeed, the coupling be- 
tween the TASEP and LK, as is was introduced above, 
produces new phenomena and extends the interest to- 
ward systems which break conservation law in a non- 
trivial way. As we shall see in the next Section, these 
features emerge already at level of properties of the mi- 
croscopic dynamics in configuration space described by 
the master equation. 

Recently a variant of our model has been suggested 
by Popkov et al. 60]. Upon supplementing the Katz- 
Lebowitz-Spohn model by Langmuir kinetics and ana- 
lyzing it within the mean-field approach similar to |15| . 
an even richer scenario for the stationary density pro- 
file is obtained that includes the emergence of localized 
downward domain walls and the appearance of several 
'shocks' separating three distinct phases. It is also noted 
in Ref. |60j | that in general it may be important to re- 
place the mean-field current by the exact current in the 
stationary state. 

In addition to its fundamental importance for non- 
equilibrium physics in general, competition between bulk 
dynamics and boundary effects are ubiquitous in nature, 
in particular biological phenomena. The TASEP has ac- 
tually been introduced in the biophysical literature as 
a model mimicking the dynamics of ribosomes moving 
along a messenger RNA chain |29j ; for generalizations of 
these studies, see the recent work in Refs. [23,|3(j. Also 
some aspects of intracellular transport show close resem- 
blance to our model. For example, processive molecular 
motors advance along cytoskeletal filaments while attach- 
ment and detachment of motors between the cytoplasm 
and the filament occur 31] . Typically kinetic rates are 
such that these motors walk a finite fraction along the 
molecular track before detaching. This falls well into 
the regime where we expect novel stationary states. Re- 
cently, it has been shown that such dynamics can be rel- 
evant for modeling the filopod growth in eukaryotic cells 
produced by motor proteins interacting within actin fil- 
aments [32. Finally, our model could also be relevant 
for studies of surface-adsorption and growth in presence 
of biased diffusion or for traffic models with bulk on-off 
ramps [33| . 

Since our paper contains a rather comprehensive dis- 
cussion of the topic, we will give a detailed outline to 
provide the reader with some guidance through the anal- 
ysis. In Section [H] we define the model by its dynamic 
rules and make a connection to its stochastic dynamics 
on a network. Though the relation between stochastic 



dynamics and networks is interesting to fully understand 
the peculiar features of the model introduced by the com- 
bination of conserved dynamics and on/off kinetics, it 
may be skipped for the first reading. We then present 
the problem in terms of a Fock space formulation and 
discuss the symmetries of the model, both key features 
for the subsequent formulation of the mean-field theory. 
In Section IIHI we briefly discuss some technical details 
of the Monte-Carlo simulation. Then follows a key sec- 
tion of the manuscript, a detailed development of the 
mean- field approximation and the resulting "Burgers" - 
like equations in the continuum limit. Here we also dis- 
cuss a series of features of these equations which will turn 
out to be crucial for the understanding of the ensuing 
density and current profiles. 

In Sect. II VI an analytic solution of the continuum equa- 
tions is derived and compared to simulation results. We 
start the discussion for the special case that on- and off- 
rates are identical. Though simpler to analyze, this case 
is somewhat artificial as it requires a fine-tuning of the 
on- and off-rates. Generically, one expects on- and off- 
rates to differ. Then the mathematical analysis becomes 
significantly more complex. We are still able to give an 
explicit analytical solution in terms of so-called Lambert 
functions, which allows us to identify a branching point 
that explains all the features of the density profiles and 
phase diagram analytically. In particular, we find a spe- 
cial point that organizes the topology of the diagram. 
In Sect. |V]we discuss the properties of the domain wall 
characterizing the phase coexistence upon changes of the 
model parameters. In particular, we show that in the 
vicinity of the special point mentioned above the domain 
wall exhibits non-analytic behavior similar to a critical 
point in continuous phase transitions. We derive the crit- 
ical exponents and the scaling related to the amplitude 
and position of the domain wall. A conclusion, Sect. 
IVII summarizes our results and provides additional ar- 
guments on the phenomenon of phase coexistence. Last, 
we discuss some discrepancies between the mean-field ap- 
proach and the simulation results and discuss a possible 
reconciliation. 



II. THE MODEL 

In this Section we are going to describe the model in 
some detail. We will also put it into the context of net- 
work theories. This will help us to pinpoint the differ- 
ences between the TASEP and LK dynamics and show 
how a model combing both aspects will lead to novel 
phenomena. Finally, we briefly review the key ideas of 
the Fock space formulation of stochastic particle dynam- 
ics. In later chapters this formulation will be used for an 
analytic discussion of the model. 
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A. Definition of the dynamic rules 

In the microscopic model we consider a finite one- 
dimensional lattice with sites labeled i — 1, ...,N (see 
Fig. 01 and lattice spacing a — L/N, where L is the total 
length of the lattice. The site i = 1 (i = N) defines the 
left (right) boundary, while the collection i = 2, ...,N — 1 
is referred to as the bulk. 

The microscopic state of the system is characterized 
by a distribution of identical particles on the lattice, i.e. 
by configurations C — {^i=i,.../v}, where each of the oc- 
cupation numbers is equal either to zero (vacancy) or 
one (particle). We impose a hard core repulsion between 
the particles, which implies that a double or higher oc- 
cupancy of sites is forbidden in the model. The full state 
space then consists of 2 N configurations. 

The statistical properties of the model are given in 
terms of the probabilities V(C,t) to find a particular con- 
figuration C = {rii} at time t. We consider the evolution 
of the probabilities V described by a master equation: 
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dV{C,t) 
dt 



J2[Wc^cV(C',t)-Wc^c>T(C,t)\. (2) 



Here, We— >c is a non-negative transition rate from con- 
figuration C to C' . As usual, master equations con- 
serve probabilities. The microscopic processes connect- 
ing two subsequent configurations are local in configura- 
tion space. Out of the possible 2 N x 2 N transitions, we 
consider only the following elementary steps connecting 
neighboring configurations: 

A) at the site i = 2, .., N — 1 a particle can jump to site 
i + 1 if unoccupied with unit rate; 

B) at the site j = la particle can enter the lattice with 
rate a if unoccupied; 

C) at the site i — N a particle can leave the lattice 
with rate (3 if occupied. 

Additionally, in the bulk we assume that a particle: 

D) can leave the lattice with site-independent detach- 
ment rate lo£>\ 

E) can fill the site (if empty) with a rate u>a by attach- 
ment. 

Processes A) to C) constitute a totally asymmetric sim- 
ple exclusion process with open boundaries El 01 B . while 
processes D) and E) define a Langmuir kinetics [14J- We 
have taken the attachment and detachment rates to be 
independent of the particle concentration in the reser- 
voir, i.e. we have assumed that the Langmuir kinetics 
on the lattice is reaction and not diffusion limited. The 
effect of diffusion in confined geometry has been studied 
in Ref . [34J . A schematic graphical representation of the 
resulting totally asymmetric exclusion model with Lang- 
muir kinetics 15] is given in Fig. El 
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exclusion process with bulk attachment and detachment 
The entrance and exit rates at the left and right end of the 
one- dimensional lattice are given by a and /3, respectively; to a 
and lod denote the local attachment and detachment rates. 



Once we know the dynamic rules of the stochastic pro- 
cess, one may introduce the notion of neighboring con- 
figurations for C and C, if they differ only by a small 
fraction (0(1/N)) of the corresponding occupation num- 
bers. This naturally leads us to a reinterpretation of the 
dynamics in terms of networks as described in the follow- 
ing subsection. 



B. Stochastic dynamics and networks 

The Markovian dynamics of the system can be repre- 
sented in terms of a network (graph), where the configu- 
rations of the stochastic process correspond to the nodes 
(vertices) of the network. Each transition allowed by the 
dynamics is represented as a directed link (edge), and 
weighted by the corresponding transition rate which can 
be read off from the dynamic rules ^4) — E). Due to the 
local dynamics the network is very dilute. A given node 
in the network is connected to a maximum number O(N) 
of nearby configurations. Nevertheless, any configuration 
can still be reached from any point within the network. 
In other words the network is connected and does not 
break into disjunct pieces. In addition, every node has 
at least one ingoing and one outgoing link. This guaran- 
tees that the system is ergodic, at least as long as N is 
finite, and all states are recurrent [oa ]. 

On such a network a distance between two different 
configurations can be defined as the minimal number of 
steps required to connect them. Note that the "architec- 
ture" of the network corresponding to a pure TASEP is 
very different from a pure LK; see Fig. 0] for an illustra- 
tion. 

The TASEP network is characterized by large fluctu- 
ations in the connectivity. Take for example the com- 
pletely filled configuration. This state can only be left 
if the particle at the right end of the lattice is ejected 
from the system. Similarly, a configuration described by 
a step function = Q(xi — x$) with a completely filled 
lattice to the left and a completely empty lattice to the 
right of xq can only be left by a single process where the 
rightmost particle is hopping forward. We call such and 
similar states "periphery states" since they are linked to 
the rest of the network by a single or only a few outgoing 
and ingoing links. This is to be contrasted with "typi- 
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TASEP network LK network 

FIG. 4: Illustration of the network architecture corresponding 
to the totally asymmetric simple exclusion process (TASEP) 
and Langmuir kinetics (LK). 

cal states" for a given density, where particles are more 
or less randomly distributed over the lattice. Then, the 
conditional probability that an empty site is in front of a 
filled site will be finite. In other words, there will be an 
extensive number of pairs (1,0) on the lattice. This im- 
plies that a typical state will be connected with an exten- 
sive number O(N) of directed ingoing and outgoing links 
to other nodes in the network. Similarly, the shortest 
path connecting two non-neighboring configurations has 
a broad length distribution. Given two randomly chosen 
sequences of occupation numbers rij = and = 1 (i.e. 
nodes) one has to ask, how many local moves of the type 
A) to C) (i.e. links) are needed to transform one sequence 
into the other. In general, there will be a distribution of 
paths connecting these nodes. The shortest connection 
may be only a few links, if local rearrangements of par- 
ticles are sufficient for matching the microscopic config- 
urations. It seems plausible that this is the case for such 
microscopic configurations, whose coarse-grained density 
profiles are identical or at least very similar. If the spa- 
tial profiles of the coarse-grained densities corresponding 
to the two configurations differ significantly, one expects 
0(N 2 ) local rearrangements to be necessary for matching 
the microscopic configurations. This is simply a conse- 
quence of particle conservation in the bulk. For example, 
to completely empty a totally filled state obviously re- 
quires 0(N 2 ) steps. In addition, distances between two 
configurations in a TASEP network can also be highly 
asymmetric. Consider a configuration C corresponding 
to a node at the periphery of the network connected to 
a configuration C Then the corresponding reverse step 
does not exist, and in order return to the configuration 
C one has to make a large loop in configuration space. 
In summary, a network corresponding to TASEP con- 
tains only directed links. A characteristic feature is its 
heterogeneity in the connectivity of nodes and distances 
between nodes. The network contains loops, many of 
which may be very long due to the conservation law in 
the bulk. 

This has to be contrasted with the architecture of a 



network corresponding to LK. Here, the connectivity of 
all nodes is independent of the particular configuration. 
Since each occupation number ni at a given site i can 
be independently changed, the number of links outgoing 
from a node is simply N. To each outgoing link there is 
an ingoing link with weights related by detailed balance. 
Moreover, any two configurations can be reached by at 
most N transitions. Since there is no conservation law, 
only local moves (particle attachment or detachment) are 
necessary. The distance of two configurations (along the 
shortest path) in a LK network is d(C, C) = J^i^ 1 \ ni ~ 
n\ | • Since the order of the necessary attachment and 
detachment processes is irrelevant the number of such 
shortest paths is highly degenerate, and depends only 
on the distance as d\. In summary, the LK network is 
not directed, very homogeneous, highly connected and 
contains many loops of all size. 

An important distinction between LK and the TASEP 
can be clearly seen if one compares the nature of the cor- 
responding stationary states. Langmuir kinetics has a 
solution described in terms of the thermodynamic equi- 
librium distribution: 

K \c\ 

nc) = WTW ^ 2 . (3) 

Here \C\ = Tii is the number of occupied sites in 

the bulk and K = loa/^d is the binding constant . Note 
that the equilibrium distribution of LK can be charac- 
terized by a Boltzmann weight upon introducing an ef- 
fective Hamiltonian TL = — fes? 1 ^^ 1 m InK. The case 
K = 1 has an interesting topological interpretation since 
the links in the LK network loose their directionality and 
the effective Hamiltonian Tt evaluates to 0. 

In contrast, the totally asymmetric exclusion process 
does not satisfy the detailed balance condition 

W c ^c V(C) = W e ^e> V{C) , 

and evolves into a non-equilibrium steady state. Actu- 
ally, if one would assume detailed balance along a closed 
directed loop in the TASEP network, one would be lead 
to the conclusion that all probabilities along the path 
have to be zero. This, in turn, would contradict the er- 
godicity of the finite system. 

The network analogy discussed above can be used to 
understand why a stochastic dynamics combining the to- 
tally asymmetric exclusion process and Langmuir kinet- 
ics is interesting and show a range of novel features not 
contained in the TASEP or LK alone. We have seen 
that the number of links necessary to connect two non- 
neighboring states in the TASEP (0(N 2 )) is much larger 
than in LK (O(N)). Then, if we take both the weights for 
hopping and the weights for attachment and detachment 
to scale the same way, LK dynamics will dominate due 
to its higher connectivity. In order to have competition, 
the weight of each LK link has to be decreased as pre- 
scribed in the Introduction such that the weighted path 
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lengths of the TASEP and LK are comparable. Yet an- 
other way to generate competition would be to only allow 
a finite (non-extensive) number of sites to cause attach- 
ment and detachment with a system size independent 
rate j3|| . The network structure of the totally asymmet- 
ric exclusion process with Langmuir kinetics also indi- 
cates why standard matrix product ansatz methods could 
be rather difficult to implement. 



C. Fock space formulation of stochastic dynamics 



averages of Eqs. Q in order to compute the time evo- 
lution of (hi(t)) one needs the corresponding averages of 
two-point correlations such as (nj_i(i)(l — hi(t))}. This 
two-point correlation function obeys itself an equation of 
motion connecting it to three-point and four-point corre- 
lation functions. Thus we are lead to an infinite hierar- 
chy of equations of motion, as is quite generally the case 
for quantum many body systems |44l |45| . To proceed 
one can then utilize standard approximation schemes of 
many body theory. 



It is sometimes convenient to formulate problems in 
stochastic particle dynamics in terms of a quantum 
Hamiltonian representation instead of a master equation. 
This formalism was developed already some time ago by 
several groups j^ll E^l ■ In the meantime it has found 
a broad range of applications (see, e.g. Ref. ^3)- We re- 
fer the reader for details to various review articles 0, ^(j 
and lecture notes 0, • 

In our case, the occupation numbers n.i(C) constitute in a 
natural way state space functions by measuring whether 
site i is occupied (rii = 1) or not (rii = 0) in configura- 
tion C. The corresponding Heisenberg equations for hi(t) 
then read 



— hi(t) = hi-i(t) l-hi(t) -hi(t) 

+ LO A l~hi(t) ~U} D fli(t) 



1 - h i+ i(t) 



(4a) 



for any site in the bulk, while for sites at the boundaries 
one obtains 



dt 



fii(t) = a 



l-ni(i) -ni(t) l-n 2 (t) 



(4b) 



—h N (t) = h N -i(t) 1 - h N (t) - f3h N {t) 
at l i 



The first line of Eq. 14af) is the usual contribution due to 
the TASEP. Introducing the current operator 

ji(t) = hi(t) 1 - n i+ i(t) , 

one can rewrite the right hand site of this line as %-\ —ji, 
which is a discrete form of the divergence of the current. 
This part defines a dynamics which satisfies particle num- 
ber conservation. The second line of Eq. I)4a|l represents 
the additional Langmuir kinetics, which acts as source 
and sink terms in the bulk. 

These equations can now be understood as equations 
of motions for a quantum many body problem. There 
are different routes to arrive at a solution. For one- 
dimensional problems there are many instances where 
exact methods are applicable Coherent state path 
integrals are useful to explore the scaling behavior at 
critical points 0, 0, One can also try to ana- 



D. Symmetries 

The system exhibits a particle-hole symmetry in the 
following sense. A jump of a particle to the right corre- 
sponds to a vacancy move by one step to the left. Sim- 
ilarly, a particle entering the system at the left bound- 
ary can be interpreted as a vacancy leaving the lattice, 
and vice versa for the right boundary. Attachment and 
detachment of particles in the bulk is mapped to detach- 
ment and attachment of vacancies, respectively. There- 
fore, one can easily verify that the transformation 



hi(t) <->• 1 - hjsr-i(t) 

CY. < — ► (3 



(5a) 
(5b) 
(5c) 



lyze the equations of motion directly |44j, |45j . By taking 



leaves Eqs. (0J invariant. Due to this property we can re- 
strict the discussion to the cases lja > and uja = , 
i.e. to K > 1 and K = 1, respectively. Eventually, for 
loa =Wd = 0, one arrives back at the TASEP respecting 
the same particle-hole symmetry described above. 



III. SIMULATIONS, MEAN-FIELD APPROXI- 
MATION AND CONTINUUM LIMIT 



In this Section we describe the Monte-Carlo simula- 
tions (MCS) and the mean-field approximation (MFA) 
we have used to compute the stationary average profile 
(fii) and the average current (ji) = (hi(l — n^+i)). 



A. Simulations 

We have performed Monte-Carlo simulations with ran- 
dom sequential updating using the dynamical rules A) — 
E) and evaluated both time and sample averages. The 
resulting profiles coincide in both averaging procedures 
for given parameters and different system sizes. In the 
simulations, stationary profiles have been obtained ei- 
ther over 10 5 time averages (with a typical time interval 
> 10 N between each step of average) or over the same 
number of samples (in the case of sample averages). 



G 



B. Mean-field approximation and continuum limit 

Averaging Eqs. Q over the stationary ensemble re- 
lates the mean occupation number to higher order corre- 
lation functions. The mean-field approximation consists 
in neglectin g th ese correlations (random phase approxi- 
mation) [4ll45l|: 

(ni(t)n i+1 (t)) = (hi{t)) (ft<+i(t)> • (6) 

Here, averages in the stationary state (•) are actually time 
independent and correspond to either sample or time av- 
erages due to the ergodicity property of the finite system. 
In this approximation the average current is given by 

(j i ) = (n i (t))(l-(n i+1 (t))). 

Once we have defined the average density at site i as 
Pi = ( n i(t)), Eq. (|4a|l results in: 

Pi-i(l- pi)- pi(l- p i+ \)+u A (l- Pi)-UDPi = 0, (7a) 

while at the boundaries, Eqs. (|4b[) . one obtains: 

«(! - Pi) - -P2) = 0, 
Pn-iQ- - Pn) - Ppn = 0. (7b) 

Note that the average density is a real number with < 
Pi < I, and Eqs. 10 form a set of N real algebraic non- 
linear relations, which can be solved numerically. 

An explicit solution of the previous equations can be 
obtained by coarse-graining the discrete lattice with lat- 
tice constant e = L/N to a continuum, i.e. considering a 
continuum limit. To simplify notation, we fix the total 
length to unity, L = l. For large systems N 3> 1, e <C 1, 
the rescaled position variable x = i/N, < x < 1, is 
quasi-continuous. An expansion of the average density 
p(x) = pi in powers of e yields: 

P (x ±e) = p(x) ± ed x p(x) + X -e 2 d 2 xP (x) + 0(e 3 ) . (8) 

Taking the scaling of the Langmuir rates, Eq. QJ, into 
account, Eqs. J7J are to leading order in e equivalent 
to the following non-linear differential equation for the 
average profile at the stationary state 

£ -dl P + (2p - l)d xP + Q A (1 -p)-Q DP = 0. (9) 

Equations (|7b|) now translate into boundary conditions 
for the density field, p(Q)=a and p(l) = 1 — 0. This can 
be interpreted as if the system at both ends is in con- 
tact with particle reservoirs of respective fixed densities 
a and 1 — (3. Note that the binding constant K remains 
unchanged in this limit. 

For finite e, the average current is written j = — ^d x p+ 
p(l — p). In the continuum limit e — ► + , this sug- 
gests that j — p(l — p) and that the current is bounded, 
j < 1/4. However, this bound holds only if the density is 
a smooth function of the position x. We shall show that 



density discontinuities can arise in the continuum limit. 
Then, for small e, these discontinuities would appear as 
rapid crossover regions where one cannot neglect the first 
order derivative term in the current definition so that the 
relation j < 1/4 needs not to be satisfied. The inequality 
can be violated also by the additional contribution aris- 
ing from current fluctuations neglected in the mean-field 
approximation; see e.g. Fig. |5|at the system boundaries. 

The equations obtained in mean-field approximation 
and the subsequent continuum limit still respect the 
particle-hole symmetry mentioned above. In terms of 
the continuous averaged density p, the symmetry now 
reads p(x) t— > 1 — p(l — x), a<->/3, il A <^fl D . Note that 
a numerical solution of the differential equation above 
necessarily uses a discretization. Using a standard al- 
gorithm for integrating differential equations, one would 
merely recover the original mean-field equations (JJJ. 

Equation has mathematical similarities to the sta- 
tionary case of a viscous Burgers equation [4(1 147L |48| 

dtp - -^ d xP + (d P j)d x p = T A - Fd- (10) 

In the Burgers equation p is identified with the fluid ve- 
locity and j is related to this velocity via j — p 2 /2. In 
our case, the hard-core interaction between particles im- 
plies a non-linear current-density relationship. As shown 
above, one finds in the continuum limit a parabolic rela- 
tion j = p(l — p). Dissipation is due to the term ed 2 p, 
while the sources represent fluxes from and to the bulk 
reservoir T A = ^(1 — p) and Tn — Q-dP- The net source 
term T A —Tu = (K + l)Cln(pi — p) is positive or negative 
depending on whether the density p is below or above the 
Langmuir isotherm, pi = K/(K + 1), expressed in terms 
of the binding constant K = 17 A j SIy). In conjunction 
with the non-linear current-density relation this implies 
that the density of the Langmuir isotherm will act like 
an "attractor" or "repellor" . If the slope of the current- 
density relation is positive, d p j > 0, and the density at 
the left end falls below the Langmuir isotherm the bulk 
reservoir will feed particles into the system. As result, 
the density grows towards p\ as one moves away from 
the boundary. In contrast, for a negative slope d p j < 0, 
i.e. for densities larger than 1/2, the density profiles are 
"repelled" from the Langmuir isotherm. The latter case 
can also be understood as an "attraction" by the Lang- 
muir isotherm if read starting from the right end of the 
system. Then, depending on whether the density at the 
right boundary is larger or smaller than pi there is a loss 
or gain of particles from the reservoir as one moves away 
from the right boundary into the bulk. This will turn 
out to be an important principle for the discussion of the 
density profiles in later sections; see e.g. Section IT V Bl 

From the analogy to fluid dynamics problems [4!j one 
expects singularities such as shocks in the density p to 
appear in the inviscid or non-dissipative limit e — > + . 
This conclusion can also be inferred by a direct inspection 
of the non-linear differential equation @ in the limit e = 
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0. It reduces to a first order differential equation, 

{2p-l)d x p+fl A (l-p)-fl D p = 0, (II) 

instead of a second order one, while the solution still has 
to satisfy two boundary conditions. Such a boundary 
value problem is apparently over-determined. However, 
we can define solutions of Eq. (|llfl respecting only one 
of the boundary conditions. Depending on whether they 
obey the boundary conditions on the left or right end 
of the lattice we call them the left solution p a and the 
right solution pp, respectively. Then, for < e <C 1 the 
full solution of Eq. 10 will be close to p a for positions 
on the left side of the system and similarly to pp on the 
right side. In general, we can not expect both solutions 
to match continuously at some point in the bulk of the 
lattice. Instead, for a large but finite system, the solu- 
tion of Eq. J3J will exhibit a rapid crossover from the left 
to the right solution. In the limit e — » + this crossover 
regime decreases in width and eventually leads to a dis- 
continuity of the average density profile at some position 
x w . Note that the discontinuity shows up only on the 
scale of the system size, i.e. in the rescaled variable x, 
whereas on the scale of the lattice spacing the crossover 
region always covers a large number of lattice sites. 

To locate the position of the discontinuity x w in the 
limit of large system sizes N ^> 1, i.e. e — > + , it is 
useful to derive a continuity equation for the current j 
and the sources J- a, J~d- Consider Eq. @ in the form 
d x j = J-a — ^d, where j = — ^d x p + p(l — p). Integrating 
over a small region of width 25x close to x w , one obtains 
j (x w + Sx) - j (x w - Sx) = $*2-lx ^ A - Fd) dx = S £ . In 
the limit e — ► + the relation simplifies to j a (x w + Sx) — 
jp(x w — dx) = 5o, where we have defined the left current 
Ja = Pa (1 - p a ) and similarly for the right current jp. 
Now, for 6x — > + , the contribution due to the sources 
iSo is of order Sx yielding the matching condition in terms 
of the left and right currents 

ja(x w ) =jp(x w ) . (12) 

The equivalent condition for the densities reads 

Pa(x W ) = 1 - Pp(x w ) . (13) 

A discontinuity of the density profile such as a domain 
wall can appear in the system depending on whether the 
previous condition is fulfilled for < x w < 1. Relation 
(|13fl . therefore, defines implicitly where a domain wall is 
located in the system. It allows to compute the domain 
wall position x w as well as its height A w = pp(x w ) — 
Pa(x w )- The domain wall separates regions of low (p < 
1/2) and high density (p > 1/2). In the ensuing phase 
diagram this will lead to an extended regime of phase 
coexistence. 

We shall see that in addition to domain walls, there 
may appear also discontinuities in the current jr33j | , which 
are located at the boundary of the system. We refer to 
them as boundary layers. 



IV. ANALYTIC SOLUTION OF THE CONTIN- 
UUM EQUATION 

In this Section, we will show in detail how one can 
treat the continuum equations, Eq. 0, analytically in 
the limit e — ► + . We shall compare these results with 
numerical solutions of Eq. © for finite e [64[ , and 
with corresponding profiles obtained from Monte-Carlo 
simulations. For the Monte-Carlo simulation the plots 
will show the average density (hi) and the average cur- 
rent (ji) = (hi(l — hi + i)). The densities and currents 
obtained from the numerical integration of the mean- 
field equations at finite e will be indicated as p £ and 
j s = —e/2d x p £ + p e (l — p s ) in the figures, respectively. 

This discussion will result in a classification of the 
possible solutions as a function of the entry and exit 
rates a and (3, the binding constant K — Qa/^d and 
the detachment rate £Id (phase diagram). Due to the 
particle-hole symmetry we can restrict ourselves to val- 
ues K > 1. Then, there are two cases to distinguish: 
K = 1 and K > 1. For K = 1 the constant density 
profile, pi = K/(K + 1), given by the Langmuir kinet- 
ics coincides with a point of particular symmetry of the 
TASEP. Indeed, for a density of p = 1/2 the system is 
dual under particle-hole exchange, the non-linear term in 
Eq. 10 vanishes, and it corresponds to a point of maxi- 
mal current |65|. It will turn out that K — 1 introduces 
particular features and requires a specific treatment and 
discussion. Since it is technically simpler we discuss this 
case first. 



A. The symmetric case: K — 1 

The mathematical analysis is simplified by the fact 
that the attachment and detachment rates are equal, 
Ha = = ^- Then Eq. Hll|) factorizes to: 

(2p-l)(9 x p-fi)=G. (14) 

The boundary conditions read p(0) = a and p(l) — l—f3. 
Note that this equation is symmetric with respect to 
particle-hole exchange. Indeed, except for the bound- 
aries, the equation is invariant under the transformation 
p(x) i— > 1 — p(l — x). This has important consequences for 
the density profiles, as will become clear in the following. 

1. The density and current profiles 

Equation (|14|) has only two basic solutions. A con- 
stant density pi{x) = 1/2 identical to the stoichiometry 
in Langmuir kinetics and also the density in the maxi- 
mal current phase of the TASEP. The other solution is 
a linear profile p = fix + C. The value of the integra- 
tion constant C depends on the boundary condition. One 
finds C a — a and Cp = 1 — /? — ft for solutions, p a (x) 
and pp (x) , matching the density at the left and the right 



8 



boundary, respectively. Depending on how the three so- 
lutions p a (x), pp{x) and pi(x) can be matched, different 
scenarios arise for the full density profile p(x). In the 
following we discuss the characteristic features of the so- 
lution in each quadrant of the a-(3 phase diagram for 
fixed fi. 

a. Lower left quadrant: a, (3 < 1/2. In this case the 
boundary conditions enforce a density less than 1/2 and 
greater that 1/2 at the left and right boundaries, respec- 
tively. This allows for a continuous density profile, where 
a constant density of p% = 1 /2 intervenes between the two 
linear solutions emerging from the left and right bound- 
aries. The corresponding positions separating the low 
density from the maximal current phase, p a (x a ) = 1/2, 
and the maximal current phase from the high density 
phase, pp(xp) = 1/2, are given by x a = (1 — 2a)/20 > 
and X/3 = (2/3 + 20 — 1) /2S7 < 1, respectively. The phase 
boundary x a — > + moves to the left upon increasing 
the entry rate a — > 1/2 and similarly xp — > 1 _ for the 
exit-rate /3 — > l/2~. Hence, depending on the values of 
the points x a and xp, one can classify the possible so- 
lutions according to the relative ordering of the phase 
boundaries: (i) x a < xp, (ii) x a — xp, and (iii) x a > xp. 

(i) The density profile is continuous and piecewisc lin- 
ear and given by 



• a 



fix 
1/2 

tl(x-l) + l-P for 



for < x < x a , 
for x a < x < xp , (15) 
xp < x < 1 . 



One observes a region of 5-phase coexistence: a low 
density phase (LD) with p{x) < 1/2 and j(x) < 1/4 
for < x < maximal current phase (MC) with 

p(x) = 1/2 and j(x) = 1/4 for x a < x < xp and high 
density phase (HD) with p{x) > 1/2 and j(x) < 1/4 for 
Xp < x < 1. For a plot of the densities and currents see 
Fig.0 

(ii) For x a = xp the width of the intermediate maximal 
current phase vanishes and the solution becomes a simple 
linear profile, continuously matching the densities of the 
LD and HD phase. 

(iii) Upon further increasing x a over xp, the interven- 
ing maximal current phase is lost and it is no longer 
possible to continuously concatenate the linear density 
profiles of the low and high density phase. There is nec- 
essarily a density discontinuity, located at a point x w 
where the currents corresponding to the right and left 
solutions match, j a (x w ) — jp{x w ). The position of the 
ensuing domain wall may be in or outside of the system. 
This leads us to further distinguish between the following 
three subcases: 

(iiii) If x w < the density profile in the bulk is 
above 1/2, i.e. in a HD phase. The profile is entirely 
described by the solution pp(x) up to a boundary layer 
at the left end. One observes that the boundary layer 
corresponds to a discontinuity in the current. The bulk 
current jp(x — > + ) does in general not match the in- 
coming particle flux a(l — a) at the left boundary (see 
Fig. EJ. 
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FIG. 5: Average density p(x) (a) and current j(x) (b) for 
parameters a = 0.4, ft = 0.4, Q. — 0.3 and K = 1. In 
this parameter range one observes a 3-phase coexistence: a 
maximal current phase is intervening between a low and high 
density phase. The profiles are computed analytically in the 
inviscid limit (dashed lines) and numerically for e = 10~ 3 
within a mean-field approximation (solid smooth line), and 
from Monte-Carlo simulations (solid wiggly line). Note that, 
within the resolution of the figures, the Monte-Carlo results 
and the numerical mean-field results can not be distinguished. 
The analytic density profile is shown for the solutions respect- 
ing the left and right boundaries conditions, p a and pp; we 
also show the Langmuir isotherm pi — 1/2. 



(iii2 ) For < x w < 1 the domain wall is within the 
system boundaries. Then the density profile connects a 
LD profile to a HD profile via a domain wall at position 
x w = (fi — a + /3)/2fi The density profile is given 
by: 



p(x) = 



fix + a for < x < x w , 

fl(x - 1) + 1 - /3 for x w < x < 1 . 



(16) 



Here we can already illustrate an important feature of 
our model. As one can infer from Fig. [71 the current 
forms a cusp at the position of the domain wall, with 
j a {x) and jp(x) being monotonically increasing and de- 
creasing functions of x, respectively. This follows directly 
from the continuum equation, Eq. (|10J) . and the density 
dependence of the source term Ta — J'D = 20(1/2 — p), 
which is positive or negative depending on whether the 
density is smaller or larger than 1/2. Hence, the domain 
wall is located at a maximum of the current. In addi- 
tion, the strict monotonicity of the current also implies 
that the domain wall is localized. A displacement of the 
domain wall to the right of x w would result in a current 
ja > jp- This in turn would increase the influx of par- 
ticles at the left boundary, which will drive the domain 
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FIG. 6: Average density p(x) (a) and current j(x) (b) for 
parameters a = 0.4, /3 = 0.1, Q = 0.3 and K = 1. We use 
the same legend as in Fig. The bulk profile is almost com- 
pletely described by the solution pp matching only the right 
boundary condition. At the left end, the bulk density does 
not match the boundary condition. As a result, a boundary 
layer appears. Only there does one find a noticeable differ- 
ence between the profiles of the Monte-Carlo simulation, the 
numerical computation at finite e and the analytic profile for 
vanishing e. 



wall back to its original position x w |67^ . 

(iii3 ) The solution for x w > 1 can be inferred by 
particle- hole symmetry from case (iiii). The low density 
profile is given by the solution p a (x) up to a boundary 
layer at the right end. 

b. Lower right quadrant: a > 1/2, (3 < 1/2. Here 
the density at both left and right boundaries is larger 
than pi = 1/2. Two different scenarios are possible. In 
the first scenario, the slope of the density profile pp(x) 
(matching the density at the right boundary) is so small 
that Pfj{x) is always larger than pi = 1/2; this requires 
fi < 1/2 - (3. Then, the bulk of the system is in the HD 
phase with a boundary layer on the left. This scenario 
is identical to the previous case (iiii), such that there is 
no qualitative change in the bulk upon crossing the line 
a = 1/2. In other words, there is no phase boundary 
and the system remains in the HD phase. In the second 
scenario, the slope fl > 1/2 — j3 such that we have a phase 
boundary between a high density and a maximal current 
phase. This solution can also be viewed as a limit of 
the 3-phase coexistence region, where for a — > l/2~ the 
phase boundary x a leaves the system through the left 
end and a boundary layer is created replacing the LD 
region (see Figs. |EIa,b)). 



FIG. 7: Average density p(x) (a) and current j(x) (b) for 
parameters a = 0.2, fi = 0.1, Q. = 0.3 and K = 1. We 
use the same legend as in Fig. Only in proximity to the 
domain wall the results from the mean-field approximation 
show deviations from the density profile obtained by Monte- 
Carlo simulation. 
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FIG. 8: (a)-(b) Average density p(x) and current j(x) for 
a — 0.8, f3 — 0.35, fl = 0.3 and K = 1. We use the same 
legend as in Fig. Except for the left boundary layer, in fig. 
(a) the analytic solution is described by the Langmuir density 
pi — 1/2 and the density pp matching the right boundary 
condition, (c)-(d) Average density p(x) and current j(x) for 
a = 0.35, f) = 0.8 and the same fl and K as before. Note that 
the curves map to those of (a)-(b) by particle- hole symmetry. 
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c. Upper left quadrant: a < 1/2, (3 > 1/2. This 
region in parameter space is obtained using particle-hole 
symmetry from the results for the lower right quadrant 
in the preceding paragraph (see Figs. |Sfc,d)). 

d. Upper right quadrant: a, (3 > 1/2. Here two 
boundary layers are formed, and the bulk of the system 
is characterized by a constant density equal to 1/2 (see 
Fig- El- This corresponds to the maximal current phase, 
which remains unchanged as compared to the TASEP 
without particle on- and off- kinetics. Note again that 
due to K — 1 the density with maximal current coin- 
cides with the Langmuir isotherm pi = 1/2. 
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0.16 . 



0.4 x 0.6 



FIG. 9: Average density p(x) (a) and current j(x) (b) for 
parameters a — 0.8, j3 = 0.8, Q — 0.3 and K = 1. We use the 
same legend as in Fig. The bulk density profile is given by 
the Langmuir density pi = 1/2 which corresponds also to the 
maximal current phase. Due to fluctuations, neglected in the 
mean-field approximation, the current profile obtained from 
the simulation exceeds the value 1/4 at each boundary. 



2. The phase diagram 

The analysis of the current and density profiles allows 
to draw cuts of the phase diagram in the (a, /3)-plane 
for fixed values of O. Note that the particle-hole symme- 
try renders all diagrams symmetric with respect to the 
diagonal a — [3. Depending on the kinetic rate O one 
can distinguish three topologies. Topologies of the phase 
diagrams change at critical values = 1/2 and 0=1; 
see Fig. HUl 

For < O < 1/2, Fig. HOf a). the phase diagram con- 
sists of seven phases. A 3-phase coexistence region LD- 
MC-HD at the center is surrounded by three 2-phase co- 
existence regions LD-HD, MC-HD and LD-MC. Pure LD, 
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FIG. 10: Cut of the phase diagram on the (a, /3)-plane in 
the mean-field approximation for K = 1 and different val- 
ues of (a) Q = 0.3, (b) Q = 0.5, (c) n = 1.0. The 
cases (a)-(c) correspond to the three different topologies of 
phase diagrams discussed in the main text. All lines repre- 
sent continuous transitions between different regions in the 
(a,(3,K = 1,£Id = const.) cut of the 4-dimensional param- 
eter space. The line parallel to the anti-diagonal is defined 
through the relation a + /3 + Q = 1. It represents the bor- 
der line where the points x a and x g (i.e. the points where 
the left and right solutions p a and pp meet the Langmuir 
isotherm pi = 1/2) coincide, x a — xp. The phase boundaries 
of the LD-HD coexistence phase, x w = and x w = 1, corre- 
spond to regions in which the domain wall is located at one of 
the system boundaries. These lines were computed by using 
the matching conditions for the currents: j a (l) = /3(1 — 13) 
and jf){Q) = a(l — a). In figure (a), we also emphasize the 
presence of the boundary layers at the left or the right end 
of the system. These are indicated with the letters "1" and 
"r", respectively. Such boundary layers remain present in the 
same regions as in figure (a) also for increasing Q. 



HD and MC phases are contiguous to the 2-phase re- 
gions. All lines between different regions represent con- 
tinuous changes in the average density p. The 3-phase 
coexistence region, and two of the 2-phase coexistence 
regions (LD-MC and MC-HD) are characterized by con- 
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tinuous density profiles. This is mainly due to the maxi- 
mal current phase with density pi(x) = 1/2. Acting as a 
" buffer" , this phase intervenes between the LD and HD 
phase or connects the LD and HD phases with the right 
and left boundary, respectively. Discontinuities only ap- 
pear as current and density discontinuities (boundary 
layers) at the system boundaries. This has to be con- 
trasted with the density profile in the coexistence region 
between the LD and HD phase. Here, a density disconti- 
nuity in the bulk (domain wall), separating both phases, 
is formed. 

It is also interesting to consider the limit f2 — > 0, as one 
expects to recover the TASEP scenario. Indeed, using the 
previous results, it is easy to show that for decreasing CI, 
the width of the 2-phase regions, as well as that of the 
3-phase region, shrinks to zero. The resulting diagram 
reproduces the well known topology of the pure TASEP 
in the mean-field approximation (5(J| . 

Upon increasing f2 up to the value 1/2, we find the first 
topology change in the phase diagram. The HD and LD 
phases gradually disappear, leaving only the 2-phase or 
3-phase coexistence regions at fl = 1/2, see Fig. llOf b). If 
f2 becomes larger than 1, the LD-HD coexistence region 
disappears; see Fig- HUT c - ). 

The Langmuir kinetics is approached for Q — > oo. 
Although the topology of the phase diagram does not 
change anymore, the phases become almost indistin- 
guishable for large kinetic rates. Here the Langmuir 
isotherm pi — 1/2 occupies most of the bulk, whereas 
the LD and HD regions are confined to a vicinity of the 
boundaries. 



B. The generic case: K > 1 

Though simpler to analyze, the previous case K = 1 is 
somewhat artificial as it requires a fine-tuning of the on- 
and off-rates. Generally, one would expect f / 1. Due 
to particle-hole symmetry we can restrict ourselves to 
K > 1. The analysis becomes significantly more complex 
since the continuum equation for the density, Eq. (|llfl . 
no longer factorizes into a simple form as for K = 1. 



1. The density and current profiles 

To proceed, it is convenient to introduce a rescaled 
density of the form 



(19) 



K + 1 , , , 
a(x) = j^-j [2p(x) - 1] - 1 



(17) 



where a = corresponds to the Langmuir isotherm pi = 
K/(K + 1). Since the density p{x) is bound within the 
interval [0,1], the rescaled density a{x) can assume values 
within the interval [-2K/(K - l),2/(K - 1)]. Then the 
continuum equation Eq. (jl 1|) simplifies to 

d x a(x) + d x In \a{x)\ = Q D ^ R + _ ^ . (18) 



Direct integrations yield 

\cr(x) \ exp(cr(x)) = Y(x) , 

where the function Y(x) is 

f (K + l) 2 
Y(x) = \ct(x )\ exp < Q D K _ 1 {% - %o) + <r(%o) 

(20) 

and <t(xq) is the value of the reduced density at the ref- 
erence point xq. In particular, the ones that match the 
boundary condition on the left or right end of the system 
are written 

Y a (x) = \a(0) | exp j il D ^+ ^ x + a(0) } (21) 



r ; i,-i : : n ( I. ) r H> {n D V^f(x-l)+a(l) 



where the boundary values cr(0) and a (I) can be written 
in terms of a and (3 using Eq. (|17f) and the boundary 
conditions p(0) = a and p(l) = 1-/3. 

Equations of the form of Eq. (|19|) appear in various 
contexts such as enzymology, population growth pro- 
cesses and hydrodynamics (see e.g. Ref. |51|)- They are 
known to have an explicit solution written in terms of a 
special function called W -function |5lj : 



a(x) 
a{x) 



W(Y{x)) , 
W(-Y{x)) 



a(x) > 
a(x) < . 



(22) 



The Lambert VF-function (see Fig. Illfl is a multi-valued 
function with two real branches, which we refer to as 
Wo(£) and W-i(£,). The branches merge at £ = — 1/e, 
where the Lambert W- function takes the value — 1. 
The first branch, Wo(£), is defined for £ > — 1/e ; 
it diverges at infinity sub-logarithmically. The second 
branch, H / _i(£), is always negative and defined in the do- 
main — 1/e < £ < 0. In the vicinity of the point £ = — 1/e 
the function W(£) behaves like a square root of £ since 
one gets d(W = W/[(l + W)g\ by the definition of the 
Lambert W-function, W{£) exp(W(£)) = £. 

Using these properties of the Lambert H^-function, the 
branch of W is selected according to the value of the 
rescaled density a. For a G [— 2K/(K — 1), —1] the rele- 
vant solution is W^\{— Y), while for a G [—1,0] one ob- 
tains W (-Y), Finally, in the interval a e [0,2/^ - 1)] 
one finds Wq(Y). 

The solutions are matched to the boundary conditions 
at the left and right ends according to the entry or exit 
rates. The left and right solutions, p a (x) and pp{x) 1 are 
then computed from the expressions in Eqs. I|23ll upon 
using the coordinate transformation given by Eq. (|17f) . 

Fig.ll2lprovides a graphical representation of the possi- 
ble set of solutions of the first order differential equation, 
Eq. (|18fl . In order to decide which one of them are ac- 
tually physically realized, one needs to go back to the 
full equation, either in its discrete form Eq. (|7a|l or its 
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branch point 




FIG. 11: The real branches Wo(0 and W-i(£) of the Lambert 
IF-function. 
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solutions are shown as thick lines in Fig. ll2f ah They are 
monotonically increasing towards the Langmuir isotherm 
pi = K/ (K + 1) > 1/2. This can be understood as a con- 
sequence of the accumulation of particles from the bulk 
reservoir via the Langmuir kinetics with increasing dis- 
tance from the left boundary. One might now expect that 
the density will finally approach the Langmuir isotherm. 
But, this is not the case. Instead, we find that the den- 
sity p a (x) never increases beyond 1/2, where the current 
reaches its largest possible value j max = 1/4. Mathemat- 
ically, this is a direct consequence of the analytic proper- 
ties of the Lambert VF-function, which has a branching 
point at a density 1/2; see Fig. II 2f a). With decreasing a 
the site where p a (x) meets the density 1/2 moves to the 
right. At a critical value of the entry rate, a c (Q,D, K), 
the branching point of the left solution p a touches the 
right boundary. 

Similarly to the discussion in the previous paragraph, 
solutions matching the right boundary condition are sta- 
ble only if (3 < 1/2. The corresponding density profiles, 
shown as thick lines in Fig. I12f b). are always in a high 
density regime, i.e. pp{x) > 1/2. If the density at the 
right boundary matches the Langmuir isotherm, the right 
solution is flat pp(%) — pi- Otherwise, the source terms 
do not cancel, leading to a net detachment/attachment 
flux such that the right density profiles decay monotoni- 
cally towards the Langmuir isotherm as one moves from 
the right boundary to the bulk. As a consequence, the 
right density pp (x) never crosses the Langmuir isotherm. 
The density profile for [3 — 1/2 is an extremal solution 
exhibiting the lowest possible density (p — 1/2) and high- 
est current (j — 1/4) at the right end, which then also 
coincides with the branching point of the Lambert W- 
function. 

In conclusion, for the left rescaled solution a a (x), an 
entry rate < a < 1/2 implies -2K/(K - 1) < a < -1. 
Hence we have according to the previous analysis: 



a a (x) = W- 1 (-Y a (x)) <0. 



(23a) 



For the right rescaled solution a {x) 1 one finds corre- 
spondingly 



cr (x) 



W o (Xp(x))>0 , 0<(3<1- Pl 
, (3 = 1- pi 



FIG. 12: Mathematical solutions for (a) the left density p a (x) 
and (b) the right density pp(x) for K = 3, Qd =0.1 and 
different values of the entry and exit rate a and The so- 
lutions which approach the Langmuir isotherm are those for 
a > P < 1/2 (thick lines). The solutions where the branch- 
ing point coincides with the right boundary are indicated by 
a c = 0.038532... and /3 C = 1/2. 



continuous version Eq. @. Analogous to the TASEP 
a solution matching the density prescribed by t he left 
boundary condition is stable only if a < 1/2 68]. Such 



W o (-Y (x))<O , I -pi <(3<l/2, 

(23b) 

where pi = K/(K + 1) is the constant density of the 
Langmuir isotherm. After the coordinate change |JTJ}, 
the general solution of the continuum mean-field equa- 
tion at e — > + , Eq. Ijlll) . is obtained by matching left 
and right solutions p a and pp. The remaining task is now 
to identify the different scenarios where domain walls and 
boundary layers appear. Such analytic results are con- 
firmed by the numerical computation at finite e. 

a. Lower left quadrant: a, (3 < 1/2. This is the only 
case where there are solutions that approach the Lang- 
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muir isotherm in the bulk and match both boundary con- 
ditions. The full density profile is obtained by finding the 
position x w where the left and right currents coincide, i.e. 
Pa{xw) = 1 — P{)(%w)- One has to consider three cases: 
(i) < x w < 1, (ii) x w < 0, and (iii) x w > 1. 

In case (i) , a domain wall is formed separating a region 
of low density on the left with a region of high density on 
the right. Depending on whether 1 — (3 is above or below 
pi, different profiles arc observed, see Figs. dj(a),(c). 
In the case (3 = 1 — pi, one obtains a flat profile of pp 
matching the value of the Langmuir isotherm We 
note again that the left and right solutions approach the 
Langmuir isotherm in the bulk. In analogy with the case 
K = 1, the domain wall is stabilized by the current pro- 
files controlled by the boundary conditions. 
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FIG. 13: Average density p(x) (a)-(c) and corresponding cur- 
rent j(x) (b)-(d) for a, f3 < | in a parameter regime showing 
phase separation. We have chosen fin = 0.1, K = 2 and (a)- 
(b) a = 0.1, P = 0.1 or (c)-(d) a = 0.3, f3 = 0.4. Solid lines 
correspond to the numerical solution of the mean-field theory 
with e = 10~ 3 . Monte-Carlo simulations are shown as solid 
wiggly line. Fhe flat dashed line represents the Langmuir 
isotherm, pi = K/(K + 1). The other dashed lines represent 
the analytic solutions given by the branches of the Lambert 
IT-functions matching the boundary conditions on the right 
and left end, respectively. For both cases (a) and (c), the 
solution matching the left boundary condition p a is given by 
the branch of the Lambert VF-function Y a ). For the 

solution matching the right boundary condition, pp, one has 
to consider the branch Wo- For (a) the branch of W has the 
argument Yp(x), while for (c) the argument is —Yp{x) (see as 
illustration also Fig. 1121 . 



In cases (ii) and (iii), one of the two phases is confined 
to the boundary. Explicitly, for (ii) the bulk is charac- 
terized by a HD with a boundary layer at the left end, 
see Fig. IMf ah Correspondingly, for (iii) the solution ex- 
hibits a LD bulk phase accompanied by a boundary layer 
on the right end side of the system, see Fig. HST a). 

b. The upper left quadrant, a < 1/2, f3 > 1/2. As 
discussed above, for (3 > 1/2 the solutions of the first 
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FIG. 14: Average density p(x) (a) and current j(x) (b) for 
a = 0.3, (3 = 0.1, n D = 0.1 and K = 2. We use the same 
legend as in Fig. 1131 Except the left boundary layer, the bulk 
density profile is given by the Lambert W-function, pp = 
W (Yp(x)). 




FIG. 15: Average density p{x) (a) and current j{x) (b) for 
a = 0.1, (3 = 0.4, Q D = 0.1 and K = 2. We use the same 
legend as in Fig. 1131 Except the right boundary layer, the 
bulk density profile is given by the Lambert W-function, p a = 
W-i(Y a (x)). 
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order differential equation, Eq. Qllfl. matching the right 
boundary condition are physically unstable. Instead, the 
actual density profile at the right boundary approaches 
the extremal solution Wo(— ^3=1/2) of the first order dif- 
ferential equation. The density difference to the bound- 
ary value is bridged by a boundary layer, which vanishes 
in the limit e — ► + . 

For the discussion of the density profiles in the up- 
per left quadrant we can simply parallel the arguments 
used for the lower left quadrant, once the right solution 
has been substituted with the extremal one. Depend- 
ing on the matching of the current, one finds again three 
cases, a LD phase, a 2-phase LD-HD coexistence and a 
HD phase. We conclude that the phases of the lower left 
quadrant extend to /3 > 1/2 with phase boundaries which 
are independent of the exit rate f3, i.e. parallel to the (3 
axis. The HD phase for (3 > 1/2 has some interesting 
features which are genuinely distinct from the HD phase 
for (3 < 1/2. The density profile in the bulk is indepen- 
dent of the entrance and exit rates, a and (3, at the left 
and right boundaries; it is given by the extremal solution 
W (—Yp =1 / 2 )- The density approaches p(L) = 1/2 and 
hence the current the maximal possible value j max = 1/4 
at the right boundary. These features are reminiscent 
of the maximal current phase for the TASEP. The only 
difference seems to be that here current and density are 
spatially varying along the system while they are con- 
stant for the TASEP. The essential characteristic in both 
cases is that the behavior of the system is determined 
by the bulk and not the boundaries. One is reminded 
of similar behavior of the Meissner phase in supercon- 
ducting materials. In the ensuing phase diagram we will 
hence indicate this regime as the "Meissner" (M) phase 
to distinguish it from the HD phase with boundary dom- 
inated density profiles [6^ . Note also that the parameter 
range for the M phase is broadened as compared to the 
maximal current phase of the TASEP. 

c. The remaining quadrants, a > 1/2. At a — 1/2, 
the system is already in the high density phase where the 
bulk profile does not match the entry rate. Increasing a 
beyond the value 1/2, therefore merely affects the bound- 
ary layer at the left end. The density profile is given by 
the right solution pp for j3 < 1/2 or the extremal one for 
(3 > 1/2 as before; for an illustration compare Fig. ITTA 
For (3 > 1/2 the same conclusion apply as in the preced- 
ing paragraph, resulting in a "Meissner" phase for the 
upper right quadrant. 

Let us conclude this subsection with some additional 
comments on boundary layers. Boundary layers arise 
from a mismatch between the bulk profile and the bound- 
ary conditions. They can bend either upwards or down- 
wards depending on whether the left or right boundary 
rates are above or below the values of the bulk solution at 
the ends. For example, in the right lower quadrant of the 
HD phase, a change from a depletion to an accumulation 
layer at the left end of the system occurs at a = pp(Q) 
for (3 < 1/2. 
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FIG. 16: Average density p(x) (a) and current j(x) (b) for 
tt D = 0.1, K = 2, a = 0.75 and (3 = 0.75. We use the same 
legend as in Fig. 1181 Except for the left and right boundary 
layers, the bulk profile obtained from the analytic mean-field 
result is given by the branch pp(x) = Wo{— Y/s(x)) of the 
Lambert VL-function computed for j3 = 1/2. 



2. The phase diagram 

We discuss the topology of the phase diagram with re- 
spect to cuts in the (a, /3)-plane for different values of K 
and £Id- We first consider the situation in which is 
fixed and K increases, starting from values slightly larger 
than unity. Figure I17f a) shows the phase diagram for 
K = 1.1. A low density (LD) phase occupies the upper 
left of the plane, while a high density (HD) and a Meiss- 
ner (M) phase are located on the right. In between there 
is a 2-phase coexistence region (LD-HD). In the coexis- 
tence phase a domain wall is localized at the point x w in 
the bulk. The boundaries of the coexistence region in the 
phase diagram are determined by those parameters where 
the domain wall hits cither the entrance, = 0, or 

the exit of the system, x w = 1. For (3 > 1/2, the density 
profile only develops a boundary layer at the right end, 
but remains unchanged in the bulk. Since the domain 
wall position becomes independent of /3, the boundaries 
of the 2-phase coexistence region become parallel to the 
axis a = 0. It is important to remark that from the an- 
alytic results the left solution p a is strictly smaller than 
1/2, except for the special point C in the phase diagram 
where /c Q (l) = f3 — 1/2. We shall see in Section IV Bl 
that in the vicinity of this point the domain wall exhibits 
critical properties. 

Upon increasing K, the LD phase progressively shrinks 
to a region close to the /3-axis, while the size of the two 
other phases increases; see Figs. I17f b) and H7f cL A 
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change of topology occurs when the LD phase collapses 
on this axis which happens upon passing a critical value 
of K. This critical value depends on S1d and can be 
computed using the expressions in Eqs. ifHfl) and A 
further increase of K results in a decrease of the extension 
of the LD-HD region in the phase diagram, see Figs. IT7I 
Eventually, for very large K the average bulk density in 
the HD and M regions approaches saturation p bu ik = 1- 

Similarly, increasing Qd at fixed K, the same topology 
change occurs, as described above; see Fig. ^] However, 
we note the different limiting behaviors for fin — » + and 
Qd — *• oo. In the first case, we are considering the limit 
of the model to the TASEP for a given binding constant 
K (although K 1). The 2-phase coexistence region 
LD-HD shrinks continuously to the line a = (3. In the 
same limit, in the upper right quadrant, a,j3> 1/2, the 
M phase approaches continuously the MC phase of the 
TASEP. For a very large detachment rate Od, the right 
boundary of the LD-HD coexistence phase approaches a 
straight line at finite entry rate a that can be computed 
from the analytic solution as equal to 1 — pi. In the same 
limit, the average density in the bulk reaches asymptoti- 
cally the value pbuik = Pi of the Langmuir isotherm. 

Eventually, one observes that all phase boundaries be- 
tween the LD-HD coexistence and the HD phase, i.e. 
where the domain is pinned at x w — 0, intersect at the 
same point J\f for any value of the detachment rate JId. 
This nodal point J\f can be evaluated as a = /3 = 1 — pi = 
1/(K + 1). At this point, indeed, the average density 
p is given by the flat profile of the Langmuir isotherm 
pi = K/(K + 1) which is obviously independent of flrj. 
As a result, the domain wall does not move from x w = 
for any value of the detachment rate £Id- Interestingly, 
one remarks that both points C and J\f approach contin- 
uously the triple point of the TASEP a = /? = 1/2 in the 
simultaneous limit flu — > + and K — > 1 + . 



V. DOMAIN WALL PROPERTIES 

The knowledge of the analytic solution in the mean- 
field approximation allows for a detailed study of the 
behavior of the domain wall height and position upon 
a change of the system parameters. While the results for 
the symmetric case K = 1 are more or less trivial, novel 
properties emerge for K > 1. In this Section, we shall 
start from the description of the domain wall behavior on 
the (a, /?)-plane of the phase diagram along trajectories 
of constant entry or exit rates, respectively. 



A. Position and amplitude of the domain wall on 
the (a, p)— plane 
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FIG. 17: Cuts of the phase diagrams on the (a, /3)-plane 
obtained by the exact solution of the stationary mean-field 
equation l)llfl in the inviscid limit e = for Qd = 0.1 and 
(a) K = 1.1, (b) K = 3.0, (c) K = 6.0. The two lines, 
corresponding to regions in which the domain wall is located 
at x w — and x w = 1, are obtained by using the matching 
conditions for the currents: = P(l — P) and ,7/3(0) = 

a(l — a). In figure (a), we emphasize several features. With 
the letters "1" and "r" we indicate the presence of boundary 
layers in the average density profile, forming at the left or 
the right end of the system, respectively. In the lower left 
quadrant, the left and right boundary layers form whenever 
the domain wall exits the system on the left and right end side. 
At the phase boundary between the HD and M phases, for 
P — 1/2, a boundary layer forms at the right end. Note that 
also in the M phase pbuik > 1/2. The presence of boundary 
layers in the different phases of the (a, /3)-plane is conserved 
upon variation of the binding constant K. The filled black 
circle represents the critical point C where the domain wall 
exhibits critical behavior; see Sect. IV Bl This critical point 
exits the plane for large values of K, accompanied with a 
topological change of the phase diagram. 



Figures HOf a.b) show the dependence of the domain 
wall position, x w , and height, A^, on the entry rate a 
along lines of constant exit rate f3. As can be inferred 
from the structure of the phase diagram presented in the 



preceding Section, for a small enough exit rate /?, a do- 
main wall can form in the bulk with a finite amplitude 
even for a vanishing entry rate, a — 0. For larger j3 1 
one observes that the domain wall builds up with a finite 
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FIG. 18: Cuts of the phase diagrams as in Fig. I17l for K — 3 
and (a) Q D = 0.01, (b) Q D = 0.05, (c) tl D = 0.2. The white 
circle corresponds to a nodal point of the system J\f defined 
by the condition a = P = 1 — pi = 1/(K + 1). Every line 
x w — crosses this point for an increasing Qd- 



height on the right boundary only above some specific 
value of a. If one regards the domain wall height as a 
kind of order parameter for the coexistence phase such a 
behavior can be termed a first order transition. This has 
to be contrasted with the case (3 = 0.5, where the domain 
wall enters the system at x w = 1 with infinitesimal height 
at a critical entry rate a = a c . In the same terminology 
this would then be a second order transition. Indeed, as 
we are going to discuss in the next subsection, the do- 
main wall exhibits critical properties at this point. In 
the phase diagram (Tig. ll7f a1') the corresponding critical 
point is indicated as C . 

In all cases, upon increasing a and hence the influx of 
particles, the domain wall changes its position continu- 
ously from the right to the left end of the system. Then, 
at some value a which depends on (3, the domain wall 
leaves the system with a finite amplitude A„. 

Similar behaviors of the position and height of the do- 
main wall is found as a function of (3 for fixed values of a; 
see Fig. Iiyf c.d). Here one finds that, upon increasing (3 
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FIG. 19: (a)-(b) Domain wall position x w and height A w as 
a function of the entrance a for different values of the exit rate 
P at Q.d = 0.1 and K — 3. At the critical point a = q c and 
/3 = 1/2 a domain wall forms at the right end of the system 
with an infinitesimal height A w . The value of the "critical" 
entry rate is a c = 0.038532... and can be written explicitly by 
using the analytic solution in the mean-field approximation, 
see Eq. 1251 . (c)-(d) Domain wall position x w and height 
A w as a function of the exit rate [3 for different values of the 
entrance rate a at Q.d = 0.1 and K = 3. For a = a c and 
/3 = 1/2 a domain wall forms at the right end of the system 
with an infinitesimal height A„. For exit rates (3 > 1/2, both 
domain wall position x w and height A w become independent 
of p. Changes in the exit rate only affect the size and shape of 
the boundary layer on the right end, but not the bulk density 
profile. 



and hence reducing the out-flux of particles the domain 
wall position x w moves continuously from the left to the 
right boundary . For small a, a domain wall is formed 
at a finite position x w and f3 = 0. For larger entry rates, 
the domain wall forms at x w = with a finite amplitude 
only for finite values of the exit rate (3. As before, the 
amplitude of the domain wall A w vanishes only for the 
critical value a — a c at [3 = 1/2. Indeed, when a > a c 
and (3 > 1/2, one notes that the domain wall position 
x w remains constant upon changing /3. As we have ex- 
plained above, this corresponds to the situation where 
the bulk profile is unaffected by a change in the exit rate 
(M phase). Only the magnitude of the boundary layer 
changes with increasing /3. 



B. Critical properties of the domain wall 

In this Section, we discuss the domain wall properties 
close to the special point C where the domain wall forms 
with infinitesimal height. The analysis will make use of 
the analytic solution in the mean-field approximation. 
We show that the domain wall emerges as a consequence 
of a bifurcation phenomenon, and calculate the resulting 
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non-analytic behavior of its height and position. 

At the point C, the analytic solution of the mean-field 
equations is described by a low density profile p(x) = 
p a (x) that not only matches the boundary conditions at 
the left but also the one at the right end; see Fig. 150] This 
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FIG. 20: Average density profiles computed analytically in 
the inviscid limit expressed in terms of Lambert W-function 
(dashed lines) and numerically for e = 10 -3 within a mean- 
field approximation (solid smooth line). Parameters are: a = 
OLc (see Eq. 1251 for the analytic expression), (3 = 0.5, Qd = 
0.1 and different values of K. The profile is entirely given by 
the left solution p a for the value of the entry rate a c , defined 
by the condition p a (x = 1) = 1/2, and (3 = 1/2. Note that 
in this case, p a matches simultaneously the left and the right 
boundary conditions. 

implies that the left and right currents also match at the 
right end of the system. Interestingly, at this position 
the current j is maximal [to| . By a small change of the 
system parameters in the 2-phase LD-HD coexistence re- 
gion the domain wall forms at the right end characterized 
by a small height, provided that (3 < 1/2. 

Using the analytic solution l)23|) one can give an explicit 
expression of the critical point C as a function of flrj 
and K . The condition that the left boundary matches 
the value 1/2 at x = 1 translates to a a (x = 1) = — 1 or, 
using the properties of the Lambert PF-function, as: 



Y a (x = l) = 1/e. 



(24) 



From the expression of the function Y a , Eq. Q2ip. and 
the initial condition ct(0) = (2a- 1)(K + 1)/(K- 1) - 1, 
one computes the value of the critical entry rate 
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From the discussion of the phase diagram and the general 
properties of the domain wall, one already infers that 
< a < 1/2 for not too large values of flD and K. 
The set 



defines a two dimensional smooth manifold in parameter 
space [critical manifold). 

In order to study the critical properties close to this 
manifold w e ap ply standard methods of bifurcation the- 
ory |S3> HE LLill ■ We consider a smooth path in the region 
of parameter space close to the critical manifold defined 
above. At some point C this path will cross the critical 
manifold. The small quantities that describe the behav- 
ior of the domain wall close to the critical manifold are 
the distance from the right end side, Sx c = 1 — x w , and 
the domain wall height, A w = pp{x w ) — p a (x w ). These 
quantities will be expressed to leading order in terms of 
the small deviations from the critical point 5a — a — a c , 
8(3 = (3 — 1/2, and similarly for 8Vlo and 8K . 

The matching condition of the left and right currents, 
p a (x ) = 1 , can be rewritten in terms of re- 

duced densities a as o a (x w ) + ap(x w ) = —2. As a con- 
sequence, the solution close to the critical point writes 
as cr at p(x w ) = — 1 =p Act, where we have introduced the 
reduced domain wall height Act as another small quan- 
tity. The relation between 5x and Act can be obtained 
by expanding the equality 



leading to 



erg exp (073) = -Yp(x w ) , 



Sx ~ (Act) 5 



(27) 



(28) 



where the prefactor can be explicitly computed and de- 
pends only on the value of the system parameters at the 
critical point C. A second relation connecting Act to 
the small distances 6a, 6j3, 6K and 6FId arises from the 
definition of the Lambert VF-function |ct| exp(<r) = Y by 
taking the ratio 



exp [a/3 - ct 



Yp(x w ) 



(29) 



(a = a c (K, riu), (3 = 1/2, K, Sl D ) 



(26) 



Ya (,X W ) 

The important observation is that the right-hand side is 
independent of the domain wall position Eq. (HJ. 

Expanding Eq. I|29|l to leading order, one obtains: 

(Act) 3 ~ 60 = A a 6a+A p (6/3) 2 +A K 6K+A n 8n D , (30) 

where 60 is a distance along a generic path that ends 
on the critical manifold. We do not consider the non- 
generic case where the the critical manifold is approached 
tangentially. Then one finds power laws different from 
those presented below for the generic case. 

As before, the coefficients A can be computed explicitly 
and shown to depend only on the rates at the critical 
point C . Interestingly, the distance 8O does not exhibit 
a linear term in 8(3. This is due to the singular behavior 
of the density profile p a {x) close to the right boundary 
at the critical point C; see Fig. [201 Combining the two 
expansion, one finds the following power laws: 

Sx - 60 2/3 A w - 50 1/3 . (31) 

The validity of these exponents is confirmed numerically 
in Figs. [2] We also checked that the amplitudes in 
the expansions l|31|l coincide with those obtained by the 
numerical data. 
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FIG. 21: Double decimal logarithmic plots of the critical be- 
havior of the domain wall height, A w , and position from the 
right end side, 1 — x w . We obtained the plot numerically 
with the program Maple, release 7, using the analytic mean- 
field solution in the vicinity of the critical point C and ap- 
plying the matching condition over the left and right cur- 
rents, j a (x m ) = jp(x m ). (a) As a function of a starting 
from the point C on the critical manifold with coordinates 
a c = 0.038532..., /3 = 1/2, tt D = 0.1 and K = 3. (b) As a 
function of /3 starting from the point C on the critical man- 
ifold with coordinates a = 0.038532..., f3 c = 1/2, Q D = 0.1 
and K = 3. (c) As a function of K from the point C on 
the critical manifold with coordinates a = 0.2, /3 = 1/2, 
Q.D = 0.051443... and K c = 3. (d) As a function of tt D 
from the point C on the critical manifold with coordinates 
a = 0.2, P = 1/2, fi D , c = 0.051443... and K = 3. The value 
of Qd,c can be easily obtained from Eq. (12 1 1 and the initial 
condition <r(0). Note the different scaling regime for the exit 
rate j3. 



C. Further properties of the domain wall position 

In this Section we discuss how the position of the do- 
main wall x w moves upon changing Qo for fixed a and 
K > 1 and a set of different values for j3. In the first 
quadrant, a,/? < 1/2, and for very small values of 
the coexistence phase is confined to a narrow strip par- 
allel to the diagonal a — (3; see Fig. I18f ah It extends to 
the quadrant a < 1/2, (3 > 1/2, where boundary layers 
form. On the other hand, for very large Qd the coex- 
istence phase corresponds to the region a < 1 — pi, see 
Sect. IIV B 21 The interesting features therefore arise in 
the region of a < 1 — pi and [3 < 1/2. We consider a path 
in the phase diagram with fixed a, (3 and K and follow 
how it intersects the phase boundaries of the 2-phase co- 
existence region LD-HD as fto is increased. From Fig. 1221 
one can distinguish three cases. 

For a < (3 the system is in a LD phase for very small 
£Id- Then, at a critical value of ftp it enters the LD-HD 
region where a domain wall forms at the right boundary. 



A further increase of £Id results in a change of the domain 
wall position to the left, asymptotically reaching the left 
boundary for very large values of the detachment rate 

n D . 

For (3 < a, the system is in the HD region for small 
CId- By increasing the detachment rate, it enters the 
LD-HD region. Differently from the previous case, the 
domain wall now forms at the left boundary, it moves to 
the right up to a maximal position x m for intermediate 
values oi Qo, and finally for large £Id it moves back to 
the left boundary with the same asymptotic behavior as 
the previous case. 

For [3 = a, the system remains in the 2-phase coexis- 
tence region LD-HD for all values of the detachment rate 
£Id- One can prove, using the analytic solution 123|l . 
that the domain wall position stays finite even in the 
limit fljj — > + and is given by 



<t(1)(1 + <t(0)) 

2(a(0Mi) - 1) 



(32) 



where <r(0) and cr(l) are the usual boundary conditions 
written in terms the model parameters; see Eq. 1171 
Interestingly, the domain wall position x w obtained for 
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FIG. 22: Domain wall position x w in logarithmic scale as a 
function of f2r> at a — 0.2, Qd ~ .051443... and K = 3 and 
different values of /3. If j3 > a the domain wall builds up from 
the right boundary, while for j3 < a from the left boundary. 
For a = j3 the domain wall approaches the position x m which 
is independent of the decreasing detachment rate Qd- At 
large f2o the domain wall position x w always moves to the 
left boundary as 1/Q.d- 

a = (3 and vanishing flu does not reduce to the value 
given by the mean-field continuum approximation in a 
pure TASEP, i.e. x w = 1/2. In order to regain the 
TASEP position, the binding constant K has to approach 
the unity. Moreover, from Eq. (|32H one find that in the 
limit a = (3 — > l/2~, the position x m is a singular func- 
tion of the binding constant close to K = 1 + . 

The previous discussion corroborates the fact that the 
Langmuir Kinetics constitutes a singular perturbation of 
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the TASEP even in the limit of small rates, yielding new 
features that are generated by the competition between 
the two dynamics. 



VI. CONCLUSIONS 

We have presented a detailed study of a model for 
driven one-dimensional transport introduced recently in 
Ref. |15| . where the dynamics of the totally asymmetric 
simple exclusion process has been supplemented by Lang- 
muir kinetics. This non-conservative process introduces 
competition between boundary and bulk dynamics sug- 
gesting a new class of models for non-equilibrium trans- 
port. The model is inspired by essential properties of in- 
tracellular transport on cytoskeletal filaments driven by 
processive motor proteins |3ll l55| . These molecular en- 
gines move unidirectionally along cytoskeletal filaments, 
and simultaneously are subject to a binding/unbinding 
kinetics between the filament and the cytoplasm. The 
processivity of the motors implies low rates of detach- 
ment. Attachment rates can be easily tuned by changing 
the concentration of motors in the cytoplasm. In partic- 
ular one may obtain very low attachment rates using a 
low volume concentration of motors. 

The non-conservative dynamics proposed introduces a 
non-trivial stationary state, with features qualitatively 
different from both the totally asymmetric simple ex- 
clusion process and Langmuir kinetics. The competing 
dynamics leads to an unexpected spatial modulation of 
the average density profile in the bulk. For extended 
regions in parameter space, we find that the density pro- 
file exhibits discontinuities on the scale of the system 
size which is characteristic for phase separation. Fur- 
thermore, the coexisting phases manifest themselves by 
a domain wall that, contrary to the TASEP, is localized 
in the bulk. In contrast to previous models 0)113) the 
localization is not induced by local defects, but arises via 
a collective phenomenon based on a microscopically ho- 
mogeneous bulk dynamics. The resulting phase diagram 
is topologically distinct from the totally asymmetric ex- 
clusion process and exhibits new phases. 

An analytic solution for the density profile has been de- 
rived in the context of a mean-field approximation in the 
continuum limit. The properties of the average density 
for different kinetic rates are encoded in the peculiar fea- 
tures of the Lambert I^-function. In particular, the dis- 
covery of a branching point is a prerequisite to rational- 
izing the behavior of the solution. The analytic solution 
has allowed to trace and study in detail the properties of 
the phase diagram. We found that the cases of equal and 
different binding rates give rise to rather distinct topolo- 
gies in the phase diagram. The limiting cases for small or 
large kinetic rates have been computed analytically. We 
have identified special points of the phase diagram which 
are the analogue of the " triple point" (viz. where all three 
phase boundaries meet) of the totally asymmetric simple 
exclusion process. There, a domain wall builds up with 



infinitesimal height at the boundary, exhibiting critical 
features characterized by unusual mean-field exponents. 
Finally, we have discussed some limiting cases in which 
the properties of the totally asymmetric simple exclusion 
process in the mean-field approximation are recovered. 

Let us give some more intuitive arguments on the do- 
main wall formation and localization. In the limit of large 
system sizes, the corresponding time-dependent version 
of Eq. which governs the dynamics of the "coarse- 
grained" density p reads 

d t p+(l-2p)d x p = F A -J r D. (33) 

One can easily see that on hydrodynamic scales the 
source contribution on the right hand side are negligi- 
ble compared to the terms related to the transport dy- 
namics. On these scales, the local dynamics is essen- 
tially described by mass conservation just as in the totally 
asymmetric exclusion process. Neglecting the source con- 
tribution, one can give an implicit analytic solution of 
Eq. (|33H by standard methods of partial differential equa- 
tions |47tl4q|. From such analysis, one can infer the mech- 
anism of the formation of density discontinuities such as 
shocks on the hydrodynamic scale, which usually build 
up in finite time. An interesting feature of the domain 
wall is also its slope S w [7lJ as a function of the sys- 
tem size N. In Fig. [23 we show the slope of the domain 
wall as obtained from Monte-Carlo simulation for grow- 
ing system size N. We find the scaling law S w ~ N v 
with the exponent rj = 0.50 ± 0.05. This result is fully 
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FIG. 23: Domain wall slope S w estimated from Monte Carlo 
simulations as a function of the system size N. Simulations 
where performed for a = 0.2, j3 = 0.6, K = 3 and Q,d = 0.1. 

compatible with the scaling exponent rj = 1/2 computed 
for a pure totally asymmetric exclusion process. Note 
that the mean-field approximation provides a wrong ex- 
ponent rjMFA — 1- The softening of the slope compared 
to mean-field was recently explained in Ref. pSo . l6fj| 
on the basis of domain wall fluctuations |57|. In this 



20 



picture, the fluctuating domain wall performs a random 
walk just like in the totally asymmetric exclusion process. 
However, since on global scale there is no mass conser- 
vation due to the Langmuir kinetics, the current is space 
dependent and drives the domain wall to an equilibrium 
position corresponding to a cusp in the current profile. 
Such domain wall behavior can be rephrased in terms of 
a random walk in the presence of a confining potential 
[o8| . This picture also suggests that the exponent for the 
slope of the domain wall is exactly given by rj = 1/2. 

It would be interesting to study how the case of equal 
rates can be obtained by a limiting procedure of the case 
with slightly different rates. The change of topology in 
the phase diagrams should be contained in the analytic 
solution, however one suspects from Eq. Q21JI that an es- 



sential singularity appears. Furthermore, one would like 
to see if possible variants of our model introduce new 
features similar to what has been done for a reference 
dynamics, i.e. the totally asymmetric exclusion process 



dyn; 

m 
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